# -*- coding: utf-8 -*-
"""Maps Replication Script.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1h_pavFtLRLvOqVBgu3Rt-2sxsjN_m3P3

Replication Material for Figures 1 and 2 in the paper "Regional inequalities and political trust in a global context" by Lisa Dellmuth, Journal of European Public Policy (2023).

Version 1.0, written 2023-02-15

**When replicating this code please cite:**
Dellmuth, L. and Jonsson, E. (2023). Estimating Subnational-level Political Trust in National and International Institutions. Paper presented at the 15th Annual Conference on the Political Economy of International Organizations (PEIO) 4-6 May 2023, San Diego. Available at: https://www.peio.me/wp-content/uploads/PEIO15/PEIO15_paper_53.pdf.


This Python-file replicates two world maps with subnational trust measures. It uses GADM data (version 4.1). Available at: https://gadm.org/download_world.html. Two Geopackage files are needed. First, download the GADM data as a single database in the GeoPackage format. Second, download the other version with six separate layers (one for each level of subdivision/aggregation), as a geopackage database.

#**Connect to Google drive**
"""

# Mount drive
from google.colab import drive
drive.mount('/content/drive', force_remount = True)

"""#**Install dependencies needed to create maps**"""

!apt install python3-cartopy

!pip install cartopy==v0.19.0.post1

!pip install folium
!pip install mapclassify
!pip install fiona
!pip install shapely
!pip install pyproj
!pip install rtree

!pip install geoplot
!pip install geopandas

!pip install mapclassify

# Commented out IPython magic to ensure Python compatibility.
import fiona
import shapely
import pyproj
import mapclassify
import rtree
import folium
import numpy as np
import descartes
import geopandas as gpd
import pandas as pd
#import geoplot as gplt
import matplotlib.pyplot as plt
import matplotlib as mpl
# %matplotlib inline

import mapclassify as mc

import matplotlib.colors as colors
import numpy as np
import matplotlib.patches as mpatches
import matplotlib.patheffects as PathEffects
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib import cm
from matplotlib.colors import ListedColormap,LinearSegmentedColormap
from matplotlib.lines import Line2D
import matplotlib.patheffects as pe
import matplotlib

# Path to files
path_to_files = '/content/drive/MyDrive/Subnational Trust Database/'

"""#**Load spatial and trust data**"""

#set font for graphs
from matplotlib import rcParams
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['DejaVu Sans']

#colormap
! pip install typhon
import typhon

#get names of layers in geopackage
layers = fiona.listlayers(path_to_files + 'Spatial Data/gadm_410-levels.gpkg')
layers

# load geopackage layers
gadm0 = gpd.read_file(path_to_files + 'Spatial Data/gadm_410-levels.gpkg', layer='ADM_0')
gadm1 = gpd.read_file(path_to_files + 'Spatial Data/gadm_410-levels.gpkg', layer='ADM_1')
gadm2 = gpd.read_file(path_to_files + 'Spatial Data/gadm_410-levels.gpkg', layer='ADM_2')
# combined geopackage for base layer
gadm = gpd.read_file(path_to_files + 'Spatial Data/gadm_410.gpkg')
#flatten map crs to see europe better
gadm = gadm.to_crs("EPSG:3395")
gadm0 = gadm0.to_crs("EPSG:3395")
gadm1 = gadm1.to_crs("EPSG:3395")
gadm2 = gadm2.to_crs("EPSG:3395")

#import subnational trust database file that has the GADM labels
data = pd.read_csv('/content/drive/MyDrive/Subnational Trust Database/STB_gadmlabels.csv')

#add all regional organisations to indicator
data.question = data.question.replace('eu_trust', 'regorg_trust')
data.question = data.question.replace('al_trust', 'regorg_trust')
data.question = data.question.replace('arableague_trust', 'regorg_trust')
data.question = data.question.replace('asean_trust', 'regorg_trust')
data.question = data.question.replace('au_trust', 'regorg_trust')
data.question = data.question.replace('bid_trust', 'regorg_trust')
data.question = data.question.replace('caf_trust', 'regorg_trust')
data.question = data.question.replace('cis_trust', 'regorg_trust')
data.question = data.question.replace('caf_trust', 'regorg_trust')
data.question = data.question.replace('gulfcoop_trust', 'regorg_trust')
data.question = data.question.replace('islcoop_trust', 'regorg_trust')
data.question = data.question.replace('mercosur_trust', 'regorg_trust')
data.question = data.question.replace('nafta_trust', 'regorg_trust')
data.question = data.question.replace('mlo_trust', 'regorg_trust')
data.question = data.question.replace('nato_trust', 'regorg_trust')
data.question = data.question.replace('sco_trust', 'regorg_trust')
data.question = data.question.replace('oas_trust', 'regorg_trust')
data.question = data.question.replace('saarc_trust', 'regorg_trust')
data.question = data.question.replace('tlc_trust', 'regorg_trust')

"""#**Create world maps**"""

#set path to export maps
path_to_files = '/content/drive/MyDrive/Subnational Trust Database/Vizualisations/Special Issue Paper Maps/'

#subset relavent data
d = data[(data["survey_code"] == "EVS") | (data["survey_code"] == "WVS")]
d = d[d['dataset'] == 'global_2var']
d = d[d['question_type'] == 'trust']
d = d[d['metric'] == 'pos'] #positivity
d = d[(d['question'] == 'regorg_trust') | (d['question'] == 'un_trust') | (d['institution'] == 'National government')] #select the relevant indicators
#get average prediction values of all 4 methods
a = d.groupby(['question','NAME_0']).mean()
a = a.reset_index()
b = d.groupby(['question','NAME_0','NAME_1']).mean()
b = b.reset_index()
c = d.groupby(['question','NAME_0','NAME_1','NAME_2']).mean()
c = c.reset_index()
#select indicator
a = a[a["question"] == "regorg_trust"]
b = b[b["question"] == "regorg_trust"]
c = c[c["question"] == "regorg_trust"]
#merge to each layer
a = a.merge(gadm0, how = 'left', left_on = 'NAME_0', right_on = 'COUNTRY')
b = b.merge(gadm1, how = 'left', on = 'NAME_1')
c = c.merge(gadm2, how = 'left', on = ['NAME_1','NAME_2'])
#transfrom to geodataframe
a = gpd.GeoDataFrame(a)
b = gpd.GeoDataFrame(b)
c = gpd.GeoDataFrame(c)
# colors for chloropleth
#cm=plt.cm.get_cmap('seismic')
cm=plt.get_cmap('difference')
cm = cm.reversed()
norm = mpl.colors.Normalize(vmin=0, vmax=1)
cm = mpl.cm.ScalarMappable(norm=norm, cmap=cm).cmap
#set boundaries
#minx, miny, maxx, maxy =  gadm.total_bounds
#west, south, east, north
minx, miny, maxx, maxy = (-20037508.34,
                          -8069073.669393,
                          20037508.34,
                          12376477.573180)
#plot figure
fig, ax = plt.subplots(1, figsize=(29, 19))
ax.axis('off')
#apply boundaries
ax.set_xlim(minx, maxx)
ax.set_ylim(miny, maxy)
# Grey base layer for missing data
base = gadm.plot(ax=ax, color = "dimgrey", legend=True, edgecolor="black", linewidth=.15,zorder = 1)
nat_borders = gadm0.plot(ax=ax, color = "none", legend = False, edgecolor="black", linewidth=.23, zorder = 9)
# Formatting
divider = make_axes_locatable(ax)
cax = divider.append_axes("bottom", pad = "3%", size="5%", aspect=.03)
# plot chloropleth layers with trust data
layer1 = a.plot(column = "value",
           ax=ax, zorder = 2,
           cmap = cm,
           norm = norm,
           alpha = 1,
           edgecolor="white",
           linewidth=.15,
           cax=cax,
           legend=False)
layer2 = b.plot(column = "value",
           ax=ax, zorder = 3,
           cmap = cm,
           norm = norm,
           alpha = 1,
           edgecolor="white",
           linewidth=.09,
           cax=cax,
           legend=False)
layer3 = c.plot(column = "value",
           ax=ax, zorder = 4,
           cmap = cm,
           norm = norm,
           alpha = 1,
           edgecolor="white",
           linewidth=.05,
           cax=cax,
           legend=False)
# Create color bar legend
cb = plt.colorbar(
    plt.cm.ScalarMappable(
        matplotlib.colors.Normalize(
            vmin = 0, vmax = 1), cmap = cm),
    orientation = 'horizontal',
    cax=cax)
cb.set_label(label = 'Percentage Confidence in Regional Organisations',size='16', labelpad=15)
cb.ax.tick_params(labelsize = '15')
cb.ax.set_xticklabels(['0%', '20%', '40%','60%','80%','100%'])
#import packages for custom legend
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib import patheffects
from matplotlib.patches import Circle, Wedge, Polygon
from matplotlib.collections import PatchCollection
import matplotlib.patches as mpatches
# Border legend
legend_elements = [mpatches.Patch(facecolor='dimgrey', edgecolor = 'black',linewidth=.15,label='No Data')]
borders_legend = ax.legend(handles=legend_elements, loc = "lower left", fontsize=15, bbox_to_anchor=(0.08, 0.01),borderpad = 1.15, fancybox=True, frameon=True)
#set title if needed
#ax.set_title('title', pad=40, fontdict={'fontsize': '40', 'fontweight' : '45'})
# save figure
fig.savefig(path_to_files + 'world_regorg_trust_pos_avg.png', dpi=500)

#subset relavent data
d = data[(data["survey_code"] == "EVS") | (data["survey_code"] == "WVS")]
d = d[d['dataset'] == 'global_2var']
d = d[d['question_type'] == 'trust']
d = d[d['metric'] == 'pos'] #positivity
d = d[(d['question'] == 'regorg_trust') | (d['question'] == 'un_trust') | (d['institution'] == 'National government')] #select the relevant indicators
#get average prediction values of all 4 methods
a = d.groupby(['question','NAME_0']).mean()
a = a.reset_index()
b = d.groupby(['question','NAME_0','NAME_1']).mean()
b = b.reset_index()
c = d.groupby(['question','NAME_0','NAME_1','NAME_2']).mean()
c = c.reset_index()
#select indicator
a = a[a["question"] == "natgov_trust"]
b = b[b["question"] == "natgov_trust"]
c = c[c["question"] == "natgov_trust"]
#merge to each layer
a = a.merge(gadm0, how = 'left', left_on = 'NAME_0', right_on = 'COUNTRY')
b = b.merge(gadm1, how = 'left', on = 'NAME_1')
c = c.merge(gadm2, how = 'left', on = ['NAME_1','NAME_2'])
#transfrom to geodataframe
a = gpd.GeoDataFrame(a)
b = gpd.GeoDataFrame(b)
c = gpd.GeoDataFrame(c)
# colors for water risk chloropleth
#cm=plt.cm.get_cmap('seismic')
cm=plt.get_cmap('difference')
cm = cm.reversed()
norm = mpl.colors.Normalize(vmin=0, vmax=1)
cm = mpl.cm.ScalarMappable(norm=norm, cmap=cm).cmap
#set boundaries
minx, miny, maxx, maxy = (-20037508.34,
                          -8069073.669393,
                          20037508.34,
                          12376477.573180)
#plot figure
fig, ax = plt.subplots(1, figsize=(29, 19))
ax.axis('off')
#apply boundaries
ax.set_xlim(minx, maxx)
ax.set_ylim(miny, maxy)
# Grey base layer for missing data
base = gadm0.plot(ax=ax, color = "dimgrey", legend=True, edgecolor="black", linewidth=.15,zorder = 1)
nat_borders = gadm0.plot(ax=ax, color = "none", legend = False, edgecolor="black", linewidth=.23, zorder = 9)
# Formatting
divider = make_axes_locatable(ax)
cax = divider.append_axes("bottom", pad = "3%", size="5%", aspect=.03)
# plot chloropleth layers with trust data
layer1 = a.plot(column = "value",
           ax=ax, zorder = 2,
           cmap = cm,
           norm = norm,
           alpha = 1,
           edgecolor="white",
           linewidth=.15,
           cax=cax,
           legend=False)
layer2 = b.plot(column = "value",
           ax=ax, zorder = 3,
           cmap = cm,
           norm = norm,
           alpha = 1,
           edgecolor="white",
           linewidth=.09,
           cax=cax,
           legend=False)
layer3 = c.plot(column = "value",
           ax=ax, zorder = 4,
           cmap = cm,
           norm = norm,
           alpha = 1,
           edgecolor="white",
           linewidth=.05,
           cax=cax,
           legend=False)
# Create color bar legend
cb = plt.colorbar(
    plt.cm.ScalarMappable(
        matplotlib.colors.Normalize(
            vmin = 0, vmax = 1), cmap = cm),
    orientation = 'horizontal',
    cax=cax)
cb.set_label(label = 'Percentage Confidence in the National Government',size='16', labelpad=15)
cb.ax.tick_params(labelsize = '15')
cb.ax.set_xticklabels(['0%', '20%', '40%','60%','80%','100%'])
#import packages for custom legend
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib import patheffects
from matplotlib.patches import Circle, Wedge, Polygon
from matplotlib.collections import PatchCollection
import matplotlib.patches as mpatches
# Border legend
legend_elements = [mpatches.Patch(facecolor='dimgrey', edgecolor = 'black',linewidth=.15,label='No Data')]
borders_legend = ax.legend(handles=legend_elements, loc = "lower left", fontsize=15, bbox_to_anchor=(0.08, 0.01),borderpad = 1.15, fancybox=True, frameon=True)
#set title if needed
#ax.set_title('title', pad=40, fontdict={'fontsize': '40', 'fontweight' : '45'})
# save figure
fig.savefig(path_to_files + 'world_natgov_trust_pos_avg.png', dpi=500)